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FOREWORD 


As it has in centuries past, the Earth’s magnetic field still plays 


a vital role in global navigation. All navigational aids or 
attitude/heading reference systems (AHRS), regardless of their 
operating principles, must speak a common language. That common 


language is in terms of the Earth’s magnetic deciinaticn. 
Consequently, magnetic-related navigational aids are integrated, 
in the form of computer hardware and software, into virtualiy every 
major weapons system of the Army, Air Force, Navy, and Marines. 
In order to maintain optimum performance, these systemsS must Le 
periodically updated with regard to the Earth’s magnetic field, 
which is a dynamic entity that changes slowly but erratically with 
time. 


For well over a century, it has been the responsibility of the U.S. 
Naval Oceanographic Office to monitor the Earth’s changing magnetic 
field and periodically report on these changes in the form of 
magnetic charts and mathematical models. For the past forty years, 
this task has involved an intensive data collection effort through 
the Project MAGNET program, which in April 1990 made the transition 
from primarily aeromagnetic surveying to satellite surveying with 
the launch of the Polar Orbiting Geomagnetic Survey (POGS) satellite. 
Follow-on satellite missions to secure data for future needs, well 
into the twenty-first century, are now being vigorously pursued. 


This report is a comprehensive summary of the cooperative effort 
between the U.S. Naval Oceanographic Office and the British 
Geological Survey in producing the 1990 Epoch World Magnetic Model, 
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SECTION 1. THE GEOMAG ALGORITHM AND THE 1990 MODEL 
1.0 Tntr ion 


The Earth’s magnetic field, as measured by a magnetic sensor on or 
above the Earth’s surface, is actually a composite of several magnetic 
fields generated by a variety of sources, which are superimposed on 
each other and which interact with each other. The most important 
of these geomagnetic sources are: 


the Earth’s fluid outer core; 
the Earth’s crust/upper mantle; 
the ionosphere; and 

the magnetosphere. 


200M 


The magnetic variation algorithm (GEOMAG) is a Fortran subroutine 
which is based on a spherical harmonic expansion of the Earth’s 
magnetic field, the coefficients of which comprise the World Magnetic 


Model (WMM). These coefficients are produced jointly by the U.S. 
Naval Oceanographic Office (NAVOCEANO)’s Geopotential Division and 
the British Geological Survey (BGS). The WMM is distributed by 


NAVOCEANO for the Defense Mapping Agency (DMA) in accordance with DMA 
Instructions 8000.1 and 8000.2. The WMMs are usually produced at 
5-year intervals and are composed of two parts: a main field model, 
which describes the Earth’s magnetic field at some base epoch, and 
a secular variation model, which accounts for the slow temporal 
variations in the main geomagnetic field from the base epoch to a 
maximum of 5 years beyond the base epoch. For example, the base epoch 
of the WMM-90 magnetic field model is 1990.0. This model is therefore 
considered valid between 1990.0 and 1995.0 and will subsequently be 
replaced at 1995.0 by the WMM-95 magnetic field model. 


It is extremely important to recognize that the WMM series of 
geomagnetic models and the charts produced from these models 
characterize only that portion of the Earth’s magnetic field which 
is generated by the Earth’s fluid outer core. The portions of the 
geomagnetic field generated by the Earth’s crust, mantle, ionosphere, 
and magnetosphere are not represented in these models. Consequently, 
a magnetic sensor such as a compass or magnetometer may observe spatial 
and temporal magnetic anomalies when referenced to the appropriate 
WMM. In particular, certain local, regional, and temporal magnetic 
declination anomalies can exceed 10 degrees. Anomalies of this 
magnitude are not common, but they do exist. Declination anomalies 
on the order of 3 or 4 degrees are not uncommon, but are of small 
spatial extent and relatively isolated. On land, spatial anomalies 
are produced by mountain ranges; ore deposits; ground which has been 
struck by lightning; geological faults; and cultural features such 
as trains, planes, tanks, railroad tracks, power lines, etc. In ocean 
areas, spatial anomalies are produced by continental margins, 
seamounts, oceanic ridges, trenches and fault zones, and ships and 
submarines. Temporal anomalies in either ocean or land areas can 
last from a few minutes to several days and are produced by icnospkeric 
and magnetospheric processes which are driven by the solar wind. 


Magnetic storms in particular can cause severe and persistent magnetic 
anomalies. Even in periods of quiet solar activity, significant 
spatial and temporal magnetic anomalies are found in the polar and 
equatorial regions of the Earth, where magnetic fields produced by 
ionospheric current systems, such as the auroral electrojets and the 
equatorial electrojet, are always present. Most of the possible 
sources of magnetic anomalies are comparatively isolated in either 
space or time. Therefore, from a global perspective, the 
root-mean-square (RMS), declination (DEC), and inclination (DIP) 
errors at sea level of the WMM are estimated to be less than 0.5 
degrees in ocean areas and less than 1.0 degrees in land areas at the 
Earth’s surface over the entire 5-year life of a particular modei. 
Also, the RMS errors at sea level of the horizontal (H) and total 
intensity (F) components of the WMM over ocean and land areas are 
estimated to be less than 200 nanoteslas (nT) over the entire 5-year 
life of the models. 


1.1 The Mathematical Model 


The Earth’s magnetic field has associated with it a geomagnetic 
potential V(r,6,0,t1), which can be expressed in spherical coordinates 
in terms of a spherical harmonic expansion of the following form: 


a+l , 
V(r,0,6,7) = ne 3(*) 2 {8am(t) COS + hyn (t) sin }P,"(8) (1) 


where the spherical coordinates (r,0,d) correspond to the radius from 
the center of the Earth, the colatitude (i.e., 90° - latitude), and 


the longitude. R, is the mean radius of the Earth; g,,(t) and h,,(t) are 
referred to as the Gauss coefficients at time tT, where tT is the time 


in years (e.g., 1987.312). P®(®) represents a particular associated 
Legendre polynomial of degree m and order m. These polynomials are 
functions of the colatitude 8. The Gauss coefficients are slowly 


varying functions of time and are expressed in the form: 


Bam(t) = 8am(T pace) + 8 am(t — Tepoce) (2a) 


han (T) ham(T epocn) + ham(® — Tepocn) ee) 


where Trpocy is the base epoch of the model, which for WMM-90 is 1990.0. 
Thus, &am(Tepocy) and h,(Tepocy) are the Gauss coefficients of the WMM at 
the model’s base epoch, while g,, and h,, (pronounced g,, dot and h,, 
dot) are tne annual rates of change of the Gauss coefficients. The 
Gauss coefficients g,,(Trocy) and &h.,,(Teeocy) and their annual rates of 
change are spherical harmonic coefricients. The Gauss coefficients 
San(Tepocn) and h,,,(Tepocy) Characterize the Earth’s main magnetic field at 


the base epoch of the model, Tgpeocy, while g,, and ham Characterize the 


secular change of the. Earth’s main magnetic field during the 5-year 
life of the model. These coefficients, up to degree and order 12 for 
the main field and up to degree and order 8 for the secular variation 
of the main field, comprise the WMM. Currently, the secular variation 
model from degree 8 through degree 12 is padded with zeros. 


The Earth’s magnetic field B(r,0,6,t) is a vector quantity having three 
components which correspond to the projection ct the magnetic field 
vector onto the three coordinate axes. Thus, B,(r,0,6,t) is that portion 
of the field pointing in the radial direction (i.e., perpendicular 
to the surface of the Earth), B,(r,6,6,t) is that portion of the field 
pointing locally due south, and B,(r,0,v,t) is that portion of the field 
pointing locally due east. The magnetic field vector can be computed 
from the geomagnetic potential by taking its gradient, thus: 


Bir,0,0,1) = -VV(r,0,0,7) (3) 


Consequently, the magnetic field components are related to the 
geomagnetic potential as follows: 


ov ,9, ’ 
B,(r,0,6,0) = - “ER09) (4a) 
1 aV(r, 6,4, 
Bir.0,0,1) = - = Meee (4b) 
1 oV(r.6,6,T 
B,(r,8,9,t) = —- rand a (4c) 


which yield the following spherical harmonic expansions: 


a+2 a 
B (r,8,o,%) = Zine 0[ | x {Ban(t) cosmo 4 h,,,(t) sinmo}P"(6) (Sa) 
Pein dP™ 
Bf{r,8,o,t1) = - 3(©) & 1Bam(t) COs > + Aag(T) sin mo } oe (5b) 
1 N R; at+2 , Z 
B,(r,8,0,t) = amo. ra 2% m{8am(t) sin md ~ h,q(t) cosmo}P, (@) (Sc) 


It must be noted that the Gauss coefficients g,,(t) and hA,,(t), as well 
as the associated Legendre polynomials and their derivatives, are 
Schmidt normalized by an international agreement (circa 1930) of the 
International Union of Geodesy and Geophysics. This particular 


normalization allows one to determine which terms of the spnerical 
harmonic model are the most significant simply by a cursory inspection 
of the model coefficients. The Schmidt-normalized associated Legendre 


Polynomials P7(®) are related to the unnormalized associated Legendre 
Polynomials P""(8) (note position of indices) by the following relation: 


P™(6) = S™P™(6) (6) 
The Schmidt normalization factors §5™ and the unnormalized associated 


Legendre Polynomials P”"(8) are computed via recurrence relationships 
as follows: 


P”(®) = 1 (7a) 
P’™(8) = sin@P*~'"-'(6) m=n#0 (7b) 
P™"(8) = cos@P*-'" -— ™pr->™ métn,n2l (7c) 
dP?(0) _ 

70 = 0 (7d) 
ap""(®) __, dp*~*"-1@) Peres e 

70 = sin@ 76 + cos6P (8) sm=n#0 (7e) 
dP™(0) _ GOO =. eg aps* 

70 = cos@ 70 sin 6P (8) eB m#n,n21 (7f£) 

where: 
— 2 2 

ew" = (n~1)°-m (8) 


and where 


ap)? 
arm (8) are 


(2n — 1)(2n -3) 


it is understood that the undefined polynomials P™'°(@) and 


to be set equal to zero. Similarly, 
1 
2n—-1 7 
{ Ss’ 1,0 n>0 
n 
4 


(9a) 


(9b) 


ae [(@-m+iV ae is for a 
= ntm are J=1 for m>1 tae 


Also, computed via recursion relations are the longitudinally dependent 
functions cos(mo) and sin(md), which are computed as follows: 


sin(imod) = O sm =0 (10a) 
cos(mo) = 1 sm =0 (10k} 
sin(m>) = sin(o)cos(m—1)6 + cos()sin(m —1)> sm >0 (10c) 
cos(mp) = cos(d)cos(m~-1)> ~ sin(d)sin(m ~1)o ym >0 (10d) 


1.2 Coordinate Transformations 


GEOMAG is intended to compute various components of the geomagnetic 
field in a geodetic coordinate system that uses the WGS-84 ellipsoid 
as the reference ellipsoid. However, the mathematical analysis in 
the previous section is based on spherical coordinates. Consequently, 
some coordinate transformations are necessary. A three-step procedure 
is required. 


a. Convert the geodetic latitude, longitude, and altitude (A,6,h) 
to spherical coordinates (r,0,). 


b. Compute the magnetic field components B,(r,6,6,t), Bj(r,9,0,t), and 
B,(r,9,9,7%). 


c. Rotate the magnetic field components from spherical 
coordinates to geodetic coordinates yielding the magnetic field 
components B,(A,6,h,t), By(A,,h,7), and BZ(A,6,h,t), which are the projections 
of the magnetic field vector B(A,0,h, 7) onto the X-north, Y-east, and 
Z-vertically down coordinates of a local rectangular coordinate system 
defined by the tangent plane to the ellipsoid which is concentric 
about the WGS-84 reference ellipsoid but which encompasses the point 


(A,o,A) . 
The transformations in step a are as follows: 


sind 


cos§ = (lla) 
VQ? cos?’ + sin? 
sin@ = V1-cos’@ 


(11b) 


where, if @ and b are respectively the semi-major and semi-minor axes 
of the WGS-84 ellipsoid: 


ti(a? = b2\ ine 2 
Ox hNa*-(a*-b*)sin"A + a (12) 
hNa?-(a?—b*)sin?h + b? 


Furthermore: 


ies 4_ 4 - 2 
r= h’+2hnVaX\a?—b)sin-r2 + @ iG =O} sinh (13) 


a*—(a*—b?)sin*X 


The transformation in step c depends on the rotation angle @ through 
which the magnetic field vector must be rotated in going from spherical 
to geodetic coordinates. This rotation angle is defined by the 
following rotations: 


cosa = {h+Va2cost*A+b*sin’X Vr (14a) 

sina = (a?—b2)cosAsind/{r Va2cos?A+b7sin’rA } (14b) 
Tt 

a = ha +8 (14¢) 


Consequently, the components of the magnetic field vector in geodetic 
coordinates may be computed as follows: 


B,(A,0,4,%) = —cosaB,(r,0,6,t)—sin a B,(r,8,6,t) (15a) 
By(A,0,h,7) = B,(r,8,6,7) (15b) 
B,(A,6,h,7) = sinaB,(r,0,o,7)—cos a B,(r,8,6, 7) (15c) 


From these rectangular components of the geomagnetic field, it is 
possible to construct all others. In particular, the following 
parameters may be computed: 


-/B2(0, 0, h, 1) + BAA, 0, A, 7) (Horizontal Intensity) (16a) 
/B2(A, 0h, ) +B2(A,0 h,v) (Total Intensity) (16b) 


B,(A,,A,7) 


B,(A,,h,%) 


(Declination) (16c) 


By (A, 9,h, 7%) meee cere 


ne me HEH ee eB eee ewe ea lClUe Ce 


BA,o h tT) = tan” \ = (Incli i n 6d 
[N\A V5, B (, h, ) ina .@) ) ( ) 
B r,0,h Tt — ‘ Grid Variat LOr 


1.3 The Computer Algorithm 


The Gauss coefficients at the base epoch, Tgpgcy, are stored in arrav 
C so that the lower half of array C is occupied by the even harmonic 
Gauss coefficients g,.(Trpocy), While the upper half of array C is 
occupied by the odd harmonic Gauss coefficients h,(Tepocy). Table 1 


illustrates the details of the storage scheme, which is equivalent 
to the following mathematical assignments: 


ax man 
Com = h m>n men 


mated 
which implies that: 


sm sn (18a) 


M1 
i>) 
5 


Sam 
Ree = Cana sm Sn,m#0 (18b) 


The annual rates of change of the Gauss coefficients are stored in 
array CD (which stands for C) so that the lower half of array CD is 


occupied by the even harmonic coefficients $,,, while the upper half 
of the array is occupied by the odd harmonic coefficients hy». Table 
2 illustrates the details of the storage scheme for array CD. It is 


essentially the same as table 1 for array C and corresponds to the 
following mathematical assignments: 


ye msn 
ez (; | (19) 
Ranst mon 
which implies that: 
Bim = Cam msn (20a) 
ae Se sm Sn,m #0 (20b) 


TABLE 1. ARRANGEMENT OF MAIN FIELD COEFFICIENTS IN ARRAY C,, 


0 1 2 3 4 5 6 7 8 9 
hs 
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10 | tmo| sis | sina | as | te | Bins | Sine | fer | Bins | fins | Ars] Ann | a 
aio ans | gna | gus | ne | gus | gine | tr | ere | tne | tre Bn | a 
to srs | gna | Bas | tne] srs | gre | tar | Arne | tne | tro] on] Ba 
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TABLE 2. ARRANGEMENT OF SECULAR VARIATION COEFFICIENTS IN ARRAY C,, 


: Pes es Pee [oe [te [i [te [| te Pe [ed 
[to | ea | ben | to | to | fs | tas | fs | fs | hos [hs | ns | fs 
jie | én | an | to | tm | tos | ha | fw | te | ho | he | fine | Fans 
jt | eu | be | to | te | te | tes | he | fn | tr [hs | nn | Har 
ftw | én | tn | to | em | to | tw | to | to | to [Fro] oh 
#re0| 8a | tina] Bias] de | dns | Bios Br | dna | Bis | Bar| | 
tro] én | @na] tna] dre drs | true fro | Bra | ns [Qo] Bf he 
[tno] ts] tra] dra] tae] ers] brs | dor | tne | Pre [ran] fron | ora 
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The numerical values of the Gauss coefficients at the base epoch and 
their corresponding annual rates of change for the WMM-90 geomagnetic 
model are listed in table 3. These numerical values are inserted 
into arrays C and CD through data statements. The base epoch of the 
moacel is also assigned through a data statement. In order to update 
the GEOMAG algorithm to a new epoch geomagnetic model such as WMM-95, 
it is necessary to replace only the data statements with the new model 
coefficients and the new base epoch. 


Important parameters in the GEOMAG routine and their mathematical 
correspondences are: 


A~a = 6378.137 km 
B~b = 6356.7523142 km 


RE~R,; = 6371.2 km 
TIME ~ % 

EPOCH ~ Trpocy 
DT ~ %t—Terpocn 
ALT ~ h 
SNORM(N,M) ~ S™ 
K(N,M) ~ *™ 
GLON ~ 


) 
SP(M) ~_ sin(md) 
CP(M) ~ cos(md) 


ST ~  sin(@) 
CT ~ cos(@) 
CA ~ cos(a) 
SA ~ sin(a) 
BR ~ 8B, 
BT ~ Be 
BP ~ B, 
BX ~ By 
BY ~ By 
BZ ~ B, 
BH ~ By 
DEC ~ Bp 
DIP ~ B, 
Tl ~ Br 


10 


TC ~ C+(t-Trocn)€ 
cD ~ C 


O20? 


Note that R, is not intended to be the mean radius of the WGS-84 
ellipsoid. It is the mean radius of a modified IAU-66 ellipsoid. 


The GEOMAG algorithm is organized into two modules, each with its own 
entry point. The first is an Initialization Module. Its purpose :5 
to compute all constants such as the recursion relation factors for 


the associated Legendre polynomials «*”, the Schmidt normalization 


factors S™, and any other parameters that do not depend on position 
Or time. The entry point for this module is GEOMAG (MAXDEG). The 
Parameter MAXDEG determines the maximum degree and order of the 
magnetic model to be used in the computations. Normally, MAXDEG = 
12, which is the maximum degree and order of the WMM series geomagnetic 
models. In order to reduce computation time, MAXDEG may be set to 
a number less than 12 (e.g., 8 or 10). However, the accuracy of the 
computed magnetic parameters is correspondingly reduced. MAXDEG must 
be set in the calling program. The second module is the Processing 
Module, which has the entry point 


GEOMG1 (ALT, GLAT, GLON, TIME, DEC, DIP, TI, GV). 


The purpose of this module is to compute the magnetic declination, 
inclination, total intensity, and grid variation of each geodetic 
position and time supplied to it. The units of the parameters in the 
argument list of the GEOMG1 entry point are as follows: 


ALT ~ kilometers (e.g., 5.314) (In) 
GLAT ~ degrees (e.g., 33.716) (In) 
GLON ~ degrees (e.g., -163.315) (In) 
TIME ~ years (e.g., 1992.427) (In) 
DEC ~ degrees (e.g., -121.734) (Out) 
DIP ~ degrees (e.g., 48.387) (Out) 
TI ~ nanoteslas (e.g., 35781.7) (Out) 
GV ~ degrees (e.g., 51.768) (Out) 


The computed magnetic field parameters are referenced to the WGS-84 
ellipsoid. The last parameter, GV, is the grid variation which is 
computed only in the polar regions (i.e., above + 55° latitude or 
below - 55° latitude). Outside of this region, a value of -999.0 is 
dummied in. It is referenced to grid north of a polar stereographic 
projection. The model is considered valid at altitudes ranging from 
sea level to 1000 km. 
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SECTION 2. THE 1990 EPOCH WORLD MAGNETIC MODEL (DERIVATION) 
2.0 verview 


There were four major data sets available for the 1933 model. These 
were: the MAGSAT satellite data collected during 1979 and 1980; the 
DE-2 satellite data collected from 1981 through 1983; Project MAGNET 
aeromagnetic data collectec between 1980 and 1990; and geomagnetic 
observatory annual magnetic means data collected between 1980 and 
1990. The global distribution of these data is illustrated in charts 
1 through 4. 


Four factors which affect the quality of the model produced and which 
influence the overall approach taken to produce the model are: 


The age of the data relative to the model epoch; 
The temporal coherence of the data; 

The spatial uniformity of the data; and 

The data density. 


aaA0 mM 


With respect to these factors, none of the four data sets are ideal. 
All four data sets, especially the satellite data sets, are dominated 
by older data. The Project MAGNET data, in addition, are neither 
temporally coherent nor spatially uniform. Furthermore, the 
observatory annual means data are sparse and suffer from severe spatial 
nonuniformity. 


The modeling objective is to create two spherical harmonic models. 
One model characterizes the Earth’s main (core-generated) magnetic 
field at the 1990.0 epoch. The other model characterizes the Earth’s 
secular (slow temporal) magnetic variations of Earth core origin for 
five years beyond the 1990 epoch. 


Given the objective and the available data, the following procedure 
was adopted: 


a. Use the observatory annual magnetic means to create two 
definitive secular variation models, the first covering the 5-year 
interval 1980 to 1985, and the second covering the 5-year interval 
1985 to 1990. These are referred to as the 1982.5 and 1987.5 definitive 
secular variation models, respectively. 


b. Use the observatory annual magnetic means to create, by 
extrapolation, one predictive secular variation model covering the 
5-year interval 1990 to 1995. It is referred to as the 1992.5 
predictive secular variation model. 


c. Use the two definitive secular variation models to push the 


satellite and aircraft magnetic field observations forward or backward, 
aS appropriate, to 1985.0. 


1§ 
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CHART 3. DE-2 DATA DISTRIBUTION 
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d. Create a 1985.0 epoch main field model using the time-adjusted 
Satellite and aircraft magnetic field observations via a weighted 
least-square inversion. 


e. Use the 1987.5 definitive secular variation m:dei to push the 
spherical harmonic coefficients of the 1985.0 epoch main field model 
forward to the 1990.0 epoch, thereby yielding the 1990.0 epoch main 
field model. 


£. Combine the 1990 epoch main field model coefficients with ti.e 
1992.5 predictive secular variation model coefficients to form the 
1990 World Magnetic Model, WMM-90. 


A by-product of this procedure is a revised 1985.0 epoch World Chart 
Model which is obtained by combining the 1985.0 main field coefficients 
generated in step d with the 1987.5 definitive secular variation model 
coefficients generated in step a. 


Zed Variation D A 1 (British Responsibility) 


The only data available for secular variation modeling are the 
observatory magnetic annual means, the first time derivative which 
provides information concerning the slow (greater than one year) rates 
of change of various components of the Earth’s main magnetic field 
at various geographic locations (roughly 200) around the world. 
Because of the sparsity and spatial nonuniformity of this data, it 
is possible to generate only a degree and order 8 spherical harmonic 
model of the secular variation. Furthermore, the predictive model 
is necessarily based on extrapolations of each magnetic component at 
each obServatory site. Examples of observatory annual means from a 
few selected sites such as Honolulu, Huancayo, Pilar, and Rude Skov, 
for the X-north, Y-east, and Z-vertically down components of the 
Earth’s magnetic field, are given in figures (la), (1b), (lc), (1d), 
(le) and (1f) through (4a), (4b), (4c), (4d), (4e) and (4f). The 
discontinuities in the field components at Honolulu are due to 
repositioning of the observatory at two separate instances. In several 
instances, the rate of change of one or more field components at an 
observatory has reversed direction over time intervals as short as 
two or three years. The sudden, unpredictable nature of the Earth’s 
field is well illustrated by these observatories. The first-order 
time derivative of these data contains magnetic field contributions 
from the Earth’s core as well as from the ionosphere and magnetosphere. 
It is difficult to remove the external field effects from these data 
because much of it is related to the solar cycle and many observatories 
do not have a sufficiently long history for a detailed analysis. 
Consequently, some external field effects are not removed from these 
data at the expense of a somewhat larger uncertainty in the secular 
variation model coefficients. 
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Three secular variation models are generated by performing a weighted, 
least-square tit of the degree and order 8 spherical harmonic model 
to the first-order time derivative of the observatory annual means. 
These models were supplied by the British Ceologica survey and are 
given in table 4. 


2.2 Main Field Data Analvsis (United States Resps....cility) 


The observatory annual magnetic means were not used in the main fieid 
modeling because those data contain, in addition to small extern. 
field contributions, some rather large local and regional magnetic 
biases of crustal origin. A detailed survey at each observatory site 
would be necessary to remove these biases. Such surveys have rarely 
been performed due to the prohibitive cost, logistics, and 
international politics involved. 


The MAGSAT data consisted of 30,473 vector magnetic field values 
selected from 401 of the first 804 orbits. Tominimize solar influences, 


the Kp magnetic index was required to be equal to, or less than, 
2(Kp$2). These orbits were individually edited by an interactive 
graphics process to delete field aligned current effects and spurious 
data. Also, the following corrections for magnetospheric effects due 
to the ring current, magnetopause Currents, and magnetotail currents 
were applied: 


B,(7,9,0,7) = - gi(t)sinO + {q)(t)cosd +5, (t) sind} cos ® (21a) 
By(r,8,0,7) = © \t)sind—s;(t)cosd (21b) 
B,(r,0,0,t) = qy(t)cos@ + {q,(t)cosd +5} (t) sind} sin® (21c) 


where the time-dependent coefficients are functions of the Disturbance 
Storm Time (Dst) index: 


q(t) = 19.69 -0.63Dst(t) (22a) 
q(t) = 0.38 -0.06Ds¢(t) (22b) 
si(t) = -2.90+0.17Dsz(t) (22c) 


These corrections are derived from the external magnetic field 
potential: 
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TABLE 4. SECULAR VARIATION MODELS (units: nanotesles/ year) 
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ViAr,9,0,7) = @ E(t J {q"(t)cosmo +5,"(t) sin m6}P7™(cos 8) (23) 
when N,,=1, via the relations: 

By = -By = : = (24a) 

By = +B, = — x (245) 

B, = -B, = v (24¢) 


A further correction takes into account the magnetic fields induced 
in the Earth by the external fields which, because of their time 
dependence and the generally low but finite conductivity of the crust 
and mantle, induce electric currents in the crust and mantle, which 


QO 
+0 


in turn generate secondary magnetic fields. These secondary fields 
are of internal origin and primarily affect the g; coefficients cf 
the internal magnetic potential: 
No on a a+ 
Viu7,9,6,71) = aL 3(2] {g,"(t) cosmo + h,"(t) sin mo }P™(cos 8) (25) 
rn=lm= 


Taking derivatives as before with n=l and m=0, the magnetic field 
corrections due to induction effects are: 


3 
B,(7,9,0,t) = —- (2 Jetesino (26a) 
B,(7,8,6,t) = 0 (26b) 
ay 
B(r,9,o,t) = - 2 [2] stercose (26c) 


where the induced part of the g) coefficient is given as: 


g(t) = 0.27q?(t) (27) 
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The external and induced magnetic field corrections given above wer 
subtracted from the MAGSAT observations. These corrections are based 
on previous analyses of MAGSAT data by Langel and Estes (1985) and 
by Quinn, Kerridge, and Barraclough (1986). 


No attempts were made to remove magnetic influences due to ionospheric 
currents such as those generated by solar quiet (SQ) currents, auroral 
electrojet currents, and equatorial electrojet currents, which are 
located below the MAGSAT orbit altitudes. These influences, theuch 
generated external to the Earth’s surface, are nevertheless part cf 
Vin(r,6,9,7) because their sources are internal to the point of 
observation. Consequently, separating core-generated fields from 
crustal and ionospherically generated fields measured by satellite 
magnetometers is difficult andis still aresearch matter. Fortunately, 
fields generated in the Earth’s crust and ionosphere are significantly 
attenuated at satellite altitudes. Therefore, errors in the main 
field model coefficiencs due to contamination of the satellite data 
by these two sources are comparatively small. 


MAGSAT was a joint National Aeronautics and Space Administration 
(NASA) /U.S. Geological Survey mission. These data were supplied by 
NASA in the form of investigator B tapes. 


The DE-2 satellite data set consisted of 5,100 data points gleaned 
from the low-altitude end of a comparatively eccentric orbit. This 
data set contained only scalar total intensity measurements of the 
Earth’s magnetic field. However, this data set exhibited a 
substantially higher rms error relative to the DGRF-80 model than the 
MAGSAT data. Consequently, for this and other reasons, the DE-2 data 
were not used in the final main field model determination. This data 
set originated with M. Sugiura of Japan (formerly of NASA). It was 
edited by J.R. Ridgeway of Science Applications Research Corporation 
and it was subsequently supplied to NAVOCEANO by Dr. Robert Langel 
of NASA. 


The Project MAGNET aeromagnetic data consisted of 338 high-level 
flights (2 15,000 feet) of vector component measurements. These data 
are routinely processed by NAVOCEANO at a 2-second sample rate and 
sent to the National Geophysical Data Center (NGDC) in Boulder, 
Colorado. A weak low-pass filter with a cut-off wavelength of 
approximately 7 km is routinely applied to this high-level data. The 
cut-off wavelength will vary slightly, depending on the average speed 
of the aircraft, which depends on prevailing wind conditions at the 
time of fiight. Typical flights last 10 to 12 hours and are flown 
at an avy age-speed of 440 km/hr. They are generally flown at night 
in orde. to minimize solar-driven external field effects which 
contribute to the Daily Variation (DV) of the Earth’s field. Project 
MAGNET flights are of long range in remote ocean areas, which precludes 
the monitoring of DV. Therefore, no explicit DV corrections are 
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made to the data. Also, the aircraft's vector magnetometer i 
calibrated at the NASA Coil Room Facility at the Goddard Space Fligh 
Center in Maryland at least once a year. 


The magnetic field observations returned from each Project MAGNET 
flight are routinely reduced in accordance with the fcellowing 
procedures: 


a. Rotate vector measurements from magnetometer coordinates 7 
instantaneous aircraft coordinates. This rotation involves only sma_.l 
misalignments relative to an imaginary coordinate system rigidly 
attached to the aircraft. 


b. Compensate in aircraft coordinates for the perturbing magnetic 
effects associated with the presence of tne aircraft by removing a 
field phenomenologically modeled as: 


- 23 ~ = ~ dBy(t 
Bc(t) = Bom + @ By(t) + B i (23) 


where the first term represents the permanent magnetic field generated 
by the remnant magnetization of the aircraft’s metal parts, the second 
term represents the field induced in the aircraft’s metal structure 


by the presence of the ambient field By(t), and the third term represents 


magnetic fields generated by eddy currents created on the aircraft’s 
metal surfaces by the aircraft’s motion through the Earth’s spatially 


varying field. Here, By(t) is the magnetic field measured by the 


magnetometer after it has been rotated into aircraft coordinates as 
indicated in step a. The compensation model contains 21 coefficients, 


3 in the vector Bp,,, 9 in the 3x3 matrix a, and 9 in the 3x3 matrix 


— 


B. 


c. Rotate the compensated field from instantaneous aircraft 
coordinates to geodetic coordinates, taking into account the 
misalignment of the inertial attitude device relative to the 
instantaneous aircraft coordinates. 


d. Visually edit the data in each flight via interactive graphics 
techniques. 


The compensation coefficients are determined by performing calibraticn 
flights at an altitude of 1,500 feet above a designated magnetic 
observatory. “These flights consist of a set of yaw, pitch, and roll 
maneuvers performed along the four cardinal headings (north, south, 
east, and west). The coefficients are then determined by a 
least-squares procedure that minimizes the squared difference between 
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the observatory field (upward continued and rotated into instantanecs us 
aircraft coordinates using the inertial attitude devices an the 
aircraft) and the field measured by the aircraft’s magnetometer. 
Using this minimization technique, the 21 compensation coefficients 
are determined simultaneously with 6 Euler angles (3 for the 
magnetometer misalignment mentioned in step a and 3 for the inertial 
attitude device misalignment mentioned in step c. The overall process 
is, therefore, nonlinear, requiring several iterations to converge. 
Note, however, that it is possible to determine only the relati-- 
misalignment between the magnetometer axes and the inertial systen 
axes. Therefore, the three magnetometer bias angles are arbitrarily 
set to zero so that, in practice, only the three inertial system bias 
angles (Euler angles) are computed. 


After compensation and editing, the aeromagnetic data were decimated 
to a 200-second sample interval (i.e., every hundredth point was 
selected), yielding 54,656 vector magnetic observations. The resulting 
Project MAGNET data set was finally converted from geodetic coordinates 
to spherical coordinates, using the coordinate transformations of the 
previous section. 


For modeling, it is desirable to have all data sets pushed forward 
or backward to a common epoch. Consequently, the MAGSAT data set, 
which was originally in spherical coordinates, was pushed forward to 
1985.0 via the 1982.5 secular variation model. The portion of the 
Project MAGNET data set collected prior to 1985.0 was pushed forward 
to 1985.0 via the 1982.5 secular variation model, while the portion 
of the Project MAGNET data set collected after 1985.0 was pushed 
backward to 1985.0 by the 1987.5 secular variation model. Subsequently, 
a revised 1985.0 epoch main field model was generated by performing 
a weighted least-squares fit of the degree and order 12 spherical 
harmonic model to the combined MAGSAT and Project MAGNET data sets. 
The resulting main field model, when combined with the 1987.5 secular 
variation model, is referred to as WC-85 (revised). These coefficients 
are listed in table 5. The 1990.0 main field model was produced by 
pushing the WC-85 (revised) main field spherical harmonic coefficients 
forward in time using the 1987.5 secular variation coefficients. The 
resulting 1990.0 epoch main field model was then combined with the 
1992.5 secular variation model to form WMM-90, the coefficients of 
which are listed in table 3 of section 1.3. 


2.3 ia) i . Mai “eld saan 


The modeling procedure used was a modification of that formulated by 
Cain et aly: (4967). The objective was to minimize the chi-square 


(y7) function 


C= M+Wet Gt Xe ee 


$0 


= mee 2 _— [2 


TABLE 5. WC-85 (REVISED) SCHMIDT NORMALIZED GAUSS COEFFICIENTS 


1 ¢ ~-29874.2 a0 aBieF 9 
bt 1 ~i904.5 5496.4 10.6 ee 
2 0 -2071.6 s) oe are) Jc 
2 1 3045.7 -2200.6 353 ae aD 
2 2 1688.7? =S0G61 .7 seas 
3 0 1294.7 ee 35 pas 
3 1 =2210:,1 -306.4 =659 ae 
2 2 1246.8 284.2 a0 =) 
3 3 832.4 #300. 7 -4.8 Spee 
4 0 9335 .0 0 0 
4 1 FBZ 60 23259 a, Be 
4 2 360.5 -247.6 -7.4 30 
4 3 -424.2 Tose ee) =| 
4 4 166.0 =29 65-5 #943 0 
a 0 S2122°5 0 8 0 
5 sf 354.0 43.7 = 4 are 
5 Z 255.2 148.7 es 0 
2 3 -94.6 -154.6 =3.3 0 
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6 0 5235 0 18 0 
6 1 63.7 =14s7 20) 0 
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6 3 -185.4 70.0 8 0 
6 4 328 -47.8 aa -1.7 
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TABLE 5. WC-85 (REVISED) SCHMIDT NORMALIZED GAUSS COEFFICIENTS (con.) 
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with respect to the 168 internal Gauss coefficients of a degree and 
order 12 spherical harmonic model, where: 


I 

” = Lw,(B,,-6,)° (30a) 
sel 

2 ‘ pd 

ke = 2 Wa(Ba — Da) (30D) 
I, 

% = Lw,(B,-d) (300) 
I, 

Xe = Dwr (Bp, — be, (30d) 


where the upper case B’s refer to the model values of their respective 
magnetic components, while the lower case b’s refer to the observed 


(measured) values of their respective magnetic components. The 
subscript i refers to a particular data point, the total number I of 
which may differ for each magnetic component. Each data point is 


weighted by a weight factor, w, which depends on several factors: 


a. Data type W,, 
MAGSAT = 1 
Project MAGNET = 1/4 


Project MAGNET observatory airswing calibrations yield rms errors on 
the order of 35 nT, while MAGSAT rms differences from degree 12 
spherical harmonic models yield rms values on the order of 9 nT. 
Consequently, the relative weight of the two data sets is taken to 


9 1 F f : . 
be =3;=;- This factor characterizes the relative quality of the two 


4 
data sets. 


b. The relative number of data points per equal area (5°x5° at 
the equator) cell; each cell was given equal weight. Therefore, data 
points corresponding to cells with more than the average number of 


points per cell, N, received less weight and vice versa. 


c. The relative rms error of data in a particular flight or orbit 


relative to the rms error, G, for all data of the corresponding data 
type (MAGSAT or Project MAGNET). 
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d. The relative rms error cf data of a specified type in ar equal 
area cell relative to all data of that type, o. 


e. The age cf the data relative to the model epo 
data collected five years away from this epech ight = 
approxtrately 1/3, While “data co; ected af. the. medec eooch get 4 
weight of 1. 


f£. Distance of geomagnetic latitude, ©, from the gecmagnet: 
equator. 


w = i for k=1,2,3(r,0,0) 
<2 ° | Oma : = é 
| O, |< 20 Mikes 0 for kA) > nel (31a) 
Wenn = O for k=1,2,3(r,8,6) 
9° . = _ 
paket ce = 1 for k=4(F) ia se 


This weighting scheme then takes the following mathematical form: 


Nin \( Sem \( Sem ) (2) 
Oem m a 
Wins = WW, = - - : } a; 
Nam /\ Semi /\ Sem ise) 


where the indices correspond to the following: 


ith - data point 

jth - equal area cell (1654 tctal) 

kth - magnetic component (r,0,0,F) 

lth - aircraft flight or satellite orbit 

mth - data type (MAGSAT, Project MAGNET) 

nth - geomagnetic latitude band (n=1 or n=2) 


The decay constant 7 was arbitrarily chosen to be 5 years, while: 


At, = % ~— Tepocn (33) 


where t, is the time of observation in years and Tgpgcy is 1985.0. 


Table 6 gives the overall rms errors of a particular magnetic component 
for each of the three separate data sets relative to the DGRF/IGRE 
series of WMMS. Table 7 lists the number of data points associated 
with each mangetic component for each sf the three data sets. Table 
8 lists the average number of data points per 5°x5° equal area ceil 
for each magnetic component for each data set. Rms statistics relative 
to the DGRF/IGRF series of models for the Project MAGNET data set are 
further broken down by Project ID and £.:yjyht number in table 9. Tue 
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TABLE 6. RMS ERRORS RELATIVE TO IGRFE/DGRF MCDELS 


MAGSAT 
Project MAGNET 


TABLE 7. NUMBER OF RECORDS 


sacse | suese | seese | sasse 


DE-2 


TABLE 8. AVERAGE NUMBER OF RECORDS PER CELL 


MAGSAT 


Project MAGNET 


ee ee 


TABLES 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 


PRO "rem 


A32-153 
A32=153 
AS32153 
AS21535 
AS2-153 
AS32153 
AB2=153 
A32-153 
A32=153 
A32-153 
ABZ-153 
C3Z2=252 
Ce 2e252 
O32=252 
Cu2-2o2 
CoeH2o2 
Ca2-252 
C3I2=252 
C32-253 
C3e=255 
632-253 
C32=293 
C32=253 
C32=253 
C32=255 
C32=253 
Cag=252 
C32-254 
C23-245 
C32=254 
C32-254 
C32-254 
C32-254 
C32-254 
C32-254 
C32-254 
C32=351 
C322351 
C32-351 
C3 2=351 
Caa= 521 
C32-351 
C3223 5) 
C32-331 
C32=351 
C3Z<3ok 


MODELS (RMS units: nT) 


JULIAN 
DAY 
129 
pw FS 
139 
143 
146 
149 
1353 
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TABLES 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 
MODELS (RMS units: nT){con.) 


JULIAN RMS RMS RMS RMS 
PROJECT ELIGHT YEAR DAY x x 2 E 

C32-351 4037 1982 332 64.3 LSe 7) 69.8 69.5 
C32-351 4038 1982 334 83.3 87.6 89.7 76.5 
C32-351 4039 1982 337 72.9 91.8 30.8 84.5 
C32-351 4040 1982 339 52.6 67.9 90.9 90.2 
C32-351 4041 1982 341 146.6 104.0 81.2 135.6 
C32-351 4042 1982 345 83.8 78.9 52.6 82.3 
C32-351 4044 1982 350 63.0 83.2 85.9 1o. 
C32-351 1098 1982 351 129.4 67.8 145.8 160.3 
C32-352 1089 1983 27 64.4 65.0 81.3 78.5 
C32-352 1090 1983 31 50.1 80.2 37.3 44.1 
C32-352 1099 1983 34 75.6 72> ° 40153 87.3 
C32-352 5045 1983 75 98.9 101.5 79.9 98.0 
C32-352 5046 1983 78 89.6 80.5 81.6 80.3 
C32-352 5047 1983 81 109.5 72.7 83.2 75.2 
032-352 5048 1983 84 212.6 106.5 145.5 130.5 
C32-352 4045 1983 87 255.4 236.2 116.9 138.3 
C32-352 4046 1983 91 39.4 76.6 56.8 45.1 
C32-352 4048 1983 96 59.2 93.9 58.7 69.6 
C32-352 1100 1983 99 31.8 35.3 Tin 64.3 
32-353 3076 1983 228 86.3 76.8 111.5 114.6 
C32-353 3077 1985 130 64.1 491.7 94.8 76.7 
C32-353 3080 1983 161 91.0 71.4 67.8 94.8 
C32-353 3081 1983 166 124.0 61.9 LASi2 136.6 
C32-353 3082 1983 169 96.1 51.1 98.2 93.4 
C32-353 3083 1983 171 59.7 121.5 76.2 67.0 
C32-353 3084 1983 175 73.2 110.1 96.5 93.4 
C32-353 3085 1983 179 91.6 100.1 103.9 114.5 
C32-353 3086 1983 183 82.7 64.7 87.7 105.9 
C32-353 3087 1983 186 44.8 100.9 64.8 58.9 
C32-451 3088 1983 325 73.4 73.7 78.4 83.8 
C32-451 3089 1983 327 92.5 66.1 60.9 96.1 
C32-451 3090 1983 330 10726 16923 -1l2.0 125.5 
C32-451 5050 1983 144 93.4 81.6 138.8 134.3 
C32-451 5051 1983 336 48.9 49.9 56.8 57.6 
C32-451 4051 1983 340 (74c0° A350. 203.2° aides 
C32-451 4052 1983 343 114.6 67.4 95.5 105.8 
C32-451 4053 1983 346 109.1 130.4 109.4 107.4 
C32-451 4054 1983 349 77.8 74.0 94.5 101.6 
C32-352 4046 1983 91 39.4 76.6 56.8 45.1 
C32-451 1103 1983 354 29-7 35.8 52.8 47.5 
C32-452 1108 1984 35 103.7 86.2 110.9 54.4 
C32-452 1106 1984 41 115.3 2123.8 124.9 122.6 
C32-452 $049 1984 45 105.9 58.8 156.2 63.5 
C32-452 3091 1984 51 53.2 265.5 242.9 56.3 
C32-452 5052 1984 64 120.5 194.8 111.9 97.0 
C32-452 4050 1984 67 168.6 175.2 80.7 64.3 
C32-452 4055 1984 73 146.1 01,9° D455 ° 101.3 
C32-452 4057 1984 76 61.8 58.1 123.0 44.3 
C32-452 4058 1984 78 86.2 68.6 127.7 55.3 
C32-452 4059 1984 84 51.3 139.6 158.6 45.2 
C32-452 4060 1984 88 28.8 138.5 98.0 38.3 
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TABLE 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO [GRF/DGRF 


PROJECT ELIGHT 


C32-452 
C32~452 
C32-453 
C32-453 
C32-453 
C32-453 
C32-453 
€32=453 
C32745935 
C32-453 
C32-453 
C327453 
C327453 
C32-453 
C32-453 
C32-453 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
C32-454 
Cazqoot 
Cs2=392 
C32-501 
Ca2-s04 
C32-951 
CSZroot 
C32=354 
C32-551 
32-552 
Cocoa 
CaZzsaoL 
C32-591 
Caer oo 
C32-55) 
Gs2s oot 
BS Pea ot 
G3esa02 
C3aZ2=o5¢ 
CAaZgH3a¢ 
O5e 7p tol 


MODELS (RMS units: nT) (con.) 


4061 
1107 
1109 
3092 
3094 
3095 
3096 
3097 
3098 
3099 
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TABLE 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 
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TABLE 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 
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TABLE 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 
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TABLE 9. PROJECT MAGNET FLIGHT STATISTICS RELATIVE TO IGRF/DGRF 
MODELSARMS units: nT) (con.) 
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tO Ene uniformity of Satellite data; the orbie-pysorooe tae eeies 
for MAGSAT and DE-2 data were taxen to be the same as for each ercire 
data set for each magnetic component. Note chat data cccusying selis 
for which there were fewer than 10 points were discarded tue to the 
presumed unreliability of their cell statistics. 
Now, the double summation expression over degree n and order m cr 
B,(r,6,0), Bo(r.6,0) and B,(r,6.6) can be converted to single summaticn 
expressions over the coefficient number / which ranges from 1 to 1** 
Then, we may write: 
168 
B,(r.8,) = 2 C0,lr.8,9) (34a) 
168 
B,(7,9, $) = »» C,Q,(r,9, >) ( 34d) 
168 
B,(7,9, ) = p> C,0,(7,8, ) (34¢) 
168 
B-(r, 9,6) = 2% CiQrilr,8,) (34d) 


where the set of coefficients {Che are the Gauss coefficients 2,, 


and h,, arbitrarily arranged so that: 


Sami W(n,m)=n(n+1)/2+m 
oe (35) 
h I(nzm)=n(n -1V/2+m+90 


am? 


This ordering then requires: 


a+? 
noo(*] cosmoP7(68) ; I(n,m)=n(n+1)/2+m 


Q,(r,9, 0) = R at2 (36a) 
non (4E] sinmoP™(6) ; [(n,m)=n(n—-1)/2+m+90 
r 

R; \" dP7(8) 

— neenr > d(n,m)=n(nt+1y2+m 

r d@ (365) 
Q,,(r,9,9) = 

R. y dP*(8) 

— ae a0 > l(n,m)=n(n - 1)/2+m+90 

r 
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R ar2 
n| Se) sin moP™(6)/sin® >; Unvm)=n(n+1)2+m 
(6.0) = 4 \) (366) 
Q,(7,9, 4 RY"? 
vn &e| cos moP7(6)/sin® [(n,m)=n(n-1y2+m+90 
r 
1 
Or, = Bar.8.9) to" 8,9)Q,,(r,8, 6) + B,(7.8,0)0,,(r, 8, o) + B,(r, 8,0)0,(7, 8.0)} (36d) 
F Vs 
These expressions are the most useful forms in which the schericai 
harmonic equations for the magnetic field components can be cast for 
a least-squares problem. 
Minimization of the chi-square function then requires that 
8 2 5 oy’ &C 
7 bal cle 
x 721 ac, ? 
be a minimum, where the symbol § means variation. This in turn 
requires: 
on? 
A = 0 = 1,...,168 38 
ac, J (38) 
Therefore, we must have: 
2 2 2 2 
Oe + xe + OX + OE te 0 jf = 1,..., 168 (39) 


a, " a, ~ a, ~ OX, 


which is a nonlinear system of 168 equations for the 168 unknown 
coefficient set {C,}}% . This system of equations is nonlinear since 


x? depends on Be which depends nonlinearly on the coefficients through 
the expression: 


B-(r,0,0) = ~VB2(r,0,0)+ Bar, 0,0) + Br, 0,0) me) 


where, B,, Bs, and B, all depend linearly on the coefficients. 


Consequently, after noting that: 


a? wa, l, 
oc p> C, py w,,Q,,(7,,8,,0,)2,,(7,,8,,9,) _ y w,.5,,0,,(7,,8;,9,) (4 pa) 


i 21 


2 168s ly 
al = 2 C, 2 We Qailr,, 9,,0,)Q0,(7,,8,,,) = 2 WabaQe(r,,8,,9,) (420) 
} = t= i= 
one 168 L 1, 
= = pas 2 Wail", 0,)24,(7,,8,,9.) 7 x Wy Dy Qe) (118,9,) (420) 
ONE 18 
Fal = = C, > we {O,(7;:,9,,0,)2,,(7,9,.0,) + Qal”,,8,,9,)00,(",8,,9,) + (Gr ss 
J = pe 
Ir 
Q,(7,,9;, 0,)Q,,(7,,9,,9,)} = yi wr,D¢,Qr, 
we have 
168 
TCO, = Rf = bl (42) 
where: 
I, ly 
Q,, =2 WO, (7581 992,710, + & We Qalr;,9;,9,)Qo,(7;,9;,9,) + (43) 


1, lp 
z We Ql", 8, 6)2,,(7;; 6,, $,) + z We {QC is 6,, 6)0,,(7,, 6, >,) + 


Qal7;,9;,9,)Qe,(7i» 8;,$;) + Q,(7,,8;,6)2,,(% 8;,;)} 
and 


I, ly 
R= © Wb, QF 8n9) + 2 Waa Qe/("8:,%) + (44) 


i, lp 
x Wybdy Qy ("is @,, 6;) + x We bp. Oeil, 6,, 6;) 
This system of 168 equations can be written in matrix form as: 


cg = Kk (45) 


which has the inverse: 


C= OK 


This is not the soluticn to the prcsblem, however, since the r: 
side of this egquaticn aiso depends on the unknown coeff ; 
That is, each element, KR, of the vector KR depends on Q,(r,6.0), which 
depends on the unknown coefficients C, in a very ncnlinear way. 


In order to solve this system of equations, we must iterate. [If p 
is the iteraticn index, then we can let: 


a aa Ge 
any = ce + 5C, EEA a) 


Then, in matrix form, we choose the following iteration scheme: 


© “lee - 0 = +3) 

Ch = O'R OS 1 on (s 
The maximum numce. TF iterations p,,, is determined by requiring that: 
168 (43) 


~S5C, < 1 nanotesla at P=Pay 
I=] 


This condition must be tested after each iteration until it is 
Satisfied. 


The rate of convergence of the algorithm depends strongly on the 
amount of noise (i.e., crustal influences, etc. in the dat 
Aeromagnetic data have a great deal of crustal noise inthem. Filtering 
the data to remove short wavelength (< 700 km) features from the data 
can improve the convergence rate by an order of magnitude. However, 
it has been shown that one-dimensional filters along the survey track 
leave short wavelength biases in the cross-track direction which 
adversely affect the final model. Consequently, no filtering (except 
for a very short wavelength (S$ 7 km) low-pass filter) was done on the 
aeromagnetic data. The number of iterations required for the model 


WAS Pmax = 14. 


In order to implement the algorithm it 1s necessary to have an initial 


guess solution C™ that is as close as possible to the actual scluticn. 
The a priori mode_ coefficients used were the existing WC-85 mode. 
eccefficients (Quinn, Kerridge, and Barraclough, 1986). 


Notice that. in. this iteration scheme, the Q matrix, which has (168 
x 168) elements (as does its inverse Q"'), needs to be computed only 
once since it does pot depend on the coefficients C;. Note, too, that 


Q iS a symmetric matrix so that only half of the elements in Q 
actually need to be computed. 
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SECTION 3... DISCUSSIGN 


4 


3.0 Meogsiine- Pesulcs 


An indication of the erratic nature of the geomagnetic field 
n ome 


jerks occur only a few times per century and are thought 
to a sudden release of magnetic energy built up fromthe electromaaqnetic 
coupling between the top of the fluid core and the lower mantle, both 
of which have substantial electrical conductivities. The numerical 
pole positions at one-year intervals for both poles are listed in 
table 10. 


Esso 

provided by the wandering of the North and South Gecmagnetic > ; 
The pole movements since 1945 are illustrated in charts 5 and 6 which 
are based on the International Geomagnetic Reference Field (IGRF) 
models, WC-85 (revised) and WMM-90. The pole movements illustrate 
a poorly understood phenomena known as the geomagnetic terk whiun 
occurred around 1970. The South Magnetic Pole movement in particulac 
illustrates a sudden change in direction at about that time These 
so ke dus 

eae 


It should be noted that roughly 90 percent of the Earth’s magnetic 
field is contained in the degree 1 spherical harmonic coefficients: 


gr, gi, and A}. These three coefficients characterize the Earth’s 


magnetic dipole field and form the basis of the geomagnetic cocrdinate 
system, which for the 1990 epoch is illustrated in chart 7. The axis 
of the geomagnetic coordinate system pierces the Earth’s surface art 
the Earth’s magnetic dipoles, which are different from the dip voles. 
The location of the dipole poles is determined when the horizontai 
(H) component of the dipole field, computed from just the degree n=] 
coefficients, is equal to zero. The dip poles on the other hand are 
determined when the horizontal (H) component of the field is computed 
using all 168 coefficients of the full-degree n=12 model is equal to 
zero. For the WMM-90 model at 1990.0, the North magnetic dipole pole 
position is located at +79.35 degrees latitude and -71.10 degrees 
longitude, while the South magnetic dipole pole position is located 
at ~-79.35 degrees latitude and +108.86 degrees longitude. The 
displacement vector for the eccentric dipole for 1990.0 in the usual 
Earth-fixed spherical coordinate system (i.e., Z-axis is the rotation 
axis, X-axis points to the prime meridian and the Y-axis is orthogonal 
to the other two, thereby creating a right-handed system) is 512 km 
radially outward from the Earth’s center, with a colatitude 21.12 
degrees and a longitude of 145.70 degrees. 


A grid of main field and annual change values of the Earth’s magnetic 
field derived from WMM-90 are tabulated in table 11 for seven basic 
magnetic field components (X,Y,Z,H,F,D,I). Contours of five of these 
components (Z,H,F,D,/) for the main field are illustrated in charts 
8 through 12. Contours of the annual change of these five components 
are illustrated in charts 13 through 17. These charts were plotted 
on a corrected Mercator projection. 
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North Pole 

Latitude Tongity 

(degrees) (degrees) 
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Contours of these same five main field magnetic components clus gris 


variation in the north polar region are given in charts 13 cnrsugs 
23, while contours of their secular variations are given in Sperce 
24 through 29. Similariy, for the south poiar region;. the main. £2324 
contours are given in charts 30 through 35). While tre. correspond: : 
secular variations are given in charts 36 through 4.1. These po- 


charts were plotted on a polar stereograpnic Progect son. . Been = 
Mercator and polar stereographic charts were generated with respec 
to the 1984 Worid Geodetic System (WGS-84) ellipsoid. 


yy bu 
ct (Db ty QQ 


Sil j nment 

The Polar Orbiting Geomagnetic Survey (POGS) satellite was launcned 
in April of 1990, too late to be used in the 1990 epoch mode:. WMM-3¢, 
having been derived from data sets independent of POGS, will Fe 4 
useful tool for evaluating POGS data and vice versa. ee ae 


quantitative comparisons between WMM-30 and the POGS data indicat 
excellent agreement between the two. The POGS data will be used = 
fabricate the 1995 epoch model. Furthermore, if the satellite remain 
operational for its maximum expected lifetime of three years, it wil 
for the first time be possible to generate a secular variation mode- 
to the same degree and order as the main field (i.e., N=M=12). 


Looking toward the end of this century and beyond, efforts have been 
made to secure data for ‘odeling purposes via the Defense Meteorological 
Satellite Program (DMSP) platform. Efforts are underway to secure 
scalar data from a boom-mounted POGS-type magnetometer on the S-:5 
DMSP satellite. This data will support the Epoch 2000 WMM. Further 
out, efforcs are being made to secure full vector magnetic capability 
from DMSP Block 6 satellites that will operate during the first quarter 


of the next century. 
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SUBROUTINE GEOMAG (GEOMAGNETIC FIELD COMPUTATION) 


eh he fe kc oie he fe he aie 2h oe ac aie ade ke ac ie fe ake 2h 2h be a ake ake ake ake ae aie afc ae akc aie ie ke kc ake 2 ee kc akc oie akc ke a i ce a 2 i fe fe ok 2k keke oi 2 kek i ai ake akc ai ke ke ode ok ake ke ke ake ok 


WMM-90 is a proposed Defense Mapping Agency (DMA) standard product. 
For information on the use and applicability of this product contact: 


DIRECTOR 

DEFENSE MAPPING AGENCY/HEADQUARTERS 
ATIN: CODE PR 

8613 LEE HIGHWAY 

FAIRFAX, VA 22031-2137 


AO ROR OR Ok OK etek rok kk kok kk kek kk ie kk ak ie ik ak ick ak kk ke ioicdc ic ac ak ack acai cick a a ak ak ak ak ioc a i akc ak cic ac ie ak ak ak akc ake akc aks 


GEOMAG PROGRAMMED BY: 


JOHN M. QUINN 7/19/90 

GEOPOTENTIAL DIVISION, CODE GGM 

U.S. NAVAL OCEANOGRAPHIC OFFICE (NAVOCEANO) 
STENNIS SPACE CENTER (SSC), MS 39522-5001 


PHONE: COM: (601) 688-4252 
AV: 485-4252 
FAX: (601) 688-5701 


Be Re Re Re He Re Re ate a RA Re ake ade he afc abe fe ote ate ake oc ec fe ake ae ae ak ak ake ae ae ae a a aie oe ake fe 2 2 ake ae ake ak ak of ak ak fe fe 2h ake 2k fe ake ade abe ofc ade abe fe fe ade ate ofe abe ate ake ke ae ate ate oe ake ake aie ak 


PURPOSE: THIS ROUTINE COMPUTES THE DECLINATION (DEC), 
INCLINATION (DIP), TOTAL INTENSITY (TD), AND 
GRID VARIATION (GV - POLAR REGIONS USING A POLAR 
STEREOGRAPHIC PROJECTION ONLY) OF THE 
EARTH’S MAGNETIC FIELD IN GEODETIC COORDINATES 
FROM THE COEFFICIENTS OF THE CURRENT OFFICIAL 
DEPARTMENT OF DEFENSE (DOD) SPHERICAL HARMONIC WORLD 
MAGNETIC MODEL (WMM-90). THE WMM SERIES OF MODELS IS 
UPDATED EVERY 5 YEARS ON | JANUARY OF THOSE YEARS 
WHICH ARE DIVISIBLE BY 5 (I.E., 1980, 1985, 1990, ETC.) 
BY THE U.S. NAVAL OCEANOGRAPHIC OFFICE IN COOPERATION 
WITH THE BRITISH GEOLOGICAL SURVEY (BGS) AND 
IS BASED ON GEOMAGNETIC SURVEY MEASUREMENTS FROM 
AIRCRAFT, SATELLITE, AND GEOMAGNETIC OBSERVATORIES. 


ANA ANANAINANANIAINIAINANIANIAANIAANNANAANANAANAANAAANAANAAANANAAANA 


eRe ae ket de teak te fe ate ate ake ate te ake de te te ake fe fe ate ate ae te ote ake ade fe ate ote oie ke oe ake i fe ee ate ake ate ke ofc ake ake ke ake 2 aie ake ok a fe fe ade fe ake ake ofc af aie ate ate 2 ie ake ate ake aie oe ade 2 ie ic fe ake ac 


MODEL: THE WMM SERIES GEOMAGNETIC MODELS ARE COMPOSED 
OF TWO PARTS: THE MAIN FIELD MODEL, WHICH IS 
VALID AT THE BASE EPOCH OF THE CURRENT MODEL, AND 
A SECULAR VARIATION MODEL, WHICH ACCOUNTS FOR SLOW 
TEMPORAL VARIATIONS IN THE MAIN GEOMAGNETIC FIELD 
FROM THE BASE EPOCH TO A MAXIMUM OF 5 YEARS BEYOND 
THE BASE EPOCH. FOR EXAMPLE, THE BASE EPOCH OF 
THE WMM-90 MODEL IS 1990.0. THIS MODEL IS THEREFORE 
CONSIDERED VALID BETWEEN 1990.0 AND 1995.0. THE 
COMPUTED MAGNETIC PARAME1ERS ARE REFERENCED TO THE 
WGS-84 ELLIPSOD. 


Bi ke ok I a RR Re a ke oi ak aie fe ac aie ke ae a ke ke ade fc 2 a ae kc oe ake ie ake a ae ioe 2c ae ke ake ake ae ac oe aie ae fe ak oie ie ake ake ake ake kc a ake ake ak ak ak a kak ok ak kk 2k ak akc akc ak a ae akc 


ACCURACY: IN OCEAN AREAS AT THE EARTH’S SURFACE OVER THE 
ENTIRE 5-YEAR LIFE OF A DEGREE AND ORDER 12 
SPHERICAL HARMONIC MODEL SUCH AS WMM-90, THE RMS 
DECLINATION ERROR IS ESTIMATED TO BE < .5 DEGREES, 
AND THE RMS INCLINATION ERROR IS ESTIMATED TO BE 
<.5 DEGREES. ALSO, THE RMS TOTAL INTENSITY ERROR 
IS ESTIMATED TO BE < 200 NANOTESLAS. 


THE ACCURACY AT ANY GIVEN TIME OF ALL FOUR 
GEOMAGNETIC PARAMETERS DEPENDS ON THE GEOMAGNETIC 
LATITUDE. THE ERRORS ARE LEAST AT THE EQUATOR AND 
GREATEST AT THE MAGNETIC POLES. 


IT IS VERY IMPORTANT TO NOTE THAT A DEGREE AND 

ORDER 12 MODEL, SUCH AS WMM-90, DESCRIBES ONLY 

THE LONG WAVELENGTH SPATIAL MAGNETIC FLUCTUATIONS 
DUE TO EARTH’S CORE. NOT INCLUDED IN THE WMM SERIES 
MODELS ARE INTERMEDIATE AND SHORT WAVELENGTH 
SPATIAL FLUCTUATIONS OF THE GEOMAGNETIC FIELD 

WHICH ORIGINATE IN THE EARTH’S MANTLE AND CRUST. 
CONSEQUENTLY, ISOLATED ANGULAR RESIDUALS AT VARIOUS 
POSITIONS ON THE SURFACE (PRIMARILY OVER LAND, IN 
CONTINENTAL MARGINS, AND OVER OCEANIC SEAMOUNTS, 
RIDGES, AND TRENCHES) OF SEVERAL DEGREES MAY BE 
EXPECTED. ALSO NOT INCLUDED IN THE MODEL ARE 
NONSECULAR TEMPORAL FLUCTUATIONS OF THE GEOMAGNETIC 
FIELD OF MAGNETOSPHERIC AND IONOSPHERIC ORIGIN. 
DURING MAGNETIC STORMS, TEMPORAL FLUCTUATIONS CAN 
CAUSE SUBSTANTIAL DEVIATIONS OF THE GEOMAGNETIC 
FIELD FROM MODEL VALUES. IN ARCTIC AND ANTARCTIC 
REGIONS, AS WELL AS IN EQUATORIAL REGIONS, DEVIATIONS 
FROM MODEL VALUES ARE BOTH FREQUENT AND PERSISTENT. 


ANNANANIAANNIANIANAINAINAINIANIANAAAAIANIAIAINAIAANAANANNAANAAANANANAANAAANN 


IF THE REQUIRED DECLINATION ACCURACY IS MORE 
STRINGENT THAN THE WMM SERIES OF MODELS PROVIDE, 
THE USER IS ADVISED TO REQUEST SPECIAL (REGIONAL OR 
LOCAL) SURVEYS BE PERFORMED AND MODELS PREPARED BY 
NAVOCEANO, WHICH OPERATES THE PROJECT MAGNET 
AIRCRAFT AND THE POLAR ORBITING GEOMAGNETIC SURVEY 
(POGS) SATELLITE. REQUESTS OF THIS NATURE SHOULD 

BE MADE THROUGH DMA AT THE ADDRESS ABOVE. 


FRO IOI IO ICI I II IK I I IORI IOI OO OR IO CI I a a aC aC aC ake aie ake fe ofc oe ofc ke ke ke Re aie ac fe kc ae ae ake ke ake 


USAGE: THIS ROUTINE IS BROKEN UP INTO TWO PARTS: 


A) AN INITIALIZATION MODULE, WHICH IS CALLED ONLY 
ONCE AT THE BEGINNING OF THE MAIN (CALLING) 
PROGRAM 

B) A PROCESSING MODULE, WHICH COMPUTES THE MAGNETIC 
FIELD PARAMETERS FOR EACH SPECIFIED GEODETIC 
POSITION (ALTITUDE, LATITUDE, LONGITUDE) AND TIME 


INITIALIZATION IS MADE VIA A SINGLE CALL TO THE MAIN 

ENTRY POINT (GEOMAG), WHILE SUBSEQUENT PROCESSING 

CALLS ARE MADE THROUGH THE SECOND ENTRY POINT (GEOMG1). 
ONE CALL TO THE PROCESSING MODULE IS REQUIRED FOR EACH 
POSITION AND TIME. 


THE VARIABLE MAXDEG IN THE INITIALIZATION CALL IS THE 
MAXIMUM DEGREE TO WHICH THE SPHERICAL HARMONIC MODEL 
IS TO BE COMPUTED. IT MUST BE SPECIFIED BY THE USER 

IN THE CALLING ROUTINE. NORMALLY IT IS 12 BUT IT MAY 

BE SET LESS THAN 12 TO INCREASE COMPUTATIONAL SPEED AT 
THE EXPENSE OF REDUCED ACCURACY. 


THE PC VERSION OF THIS SUBROUTINE MUST BE COMPILED 
WITH A FORTRAN 77 COMPATIBLE COMPILER SUCH AS THE 
MICROSOFT OPTIMIZING FORTRAN COMPILER VERSION 4.1 
OR LATER. 


FRO ROO ok igiOi kk RG GR OO Ga a ck ak ae ak ak ake ie keke ake ake ke ate ake ake oe a kk ke kak ke a oe 


REFERENCES: 


JOHN M. QUINN, DAVID J. KERRIDGE, AND DAVID R. BARRACLOUGH, 
WORLD MAGNETIC CHARTS FOR 1985 - SPHERICAL HARMONIC 
MODELS OF THE GEOMAGNETIC FIELD AND ITS SECULAR 


VARIATION, GEOPHYS. J. R. ASTR. SOC. (1986) Vol. 87, 
PP. 1143-1157 
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AAN ANIA NAAN IANINA AIAN ANN ANNAAQANNAANANAAANIAIANIAIAAAANAAAANAANQNANANANNNANAON 


DEFENSE MAPPING AGENCY TECHNICAL REPORT TR 8350.2: 
DEPARTMENT OF DEFENSE WORLD GEODETIC SYSTEM 1984, 
SEPT. 30 (1987) 


JOSEPH C. CAIN, ET AL., A PROPOSED MODEL FOR THE 
INTERNATIONAL GEOMAGNETIC REFERENCE FIELD - 1965, 
J. GEOMAG, AND GEOELECT, VOL. 19, NO. 4, PP. 335-355 
(1967) (SEE APPENDIX) 


ALFRED J. ZMUDA, WORLD MAGNETIC SURVEY 1957-1969, 
INTERNATIONAL ASSOCIATION OF GEOMAGNETISM AND 
AERONOMY (IAGA) BULLETIN #28, PP. 186-188 (1971) 
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PARAMETER DESCRIPTIONS: 
A - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM) 
B - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM) 
RE - MEAN RADIUS OF IAU-66 ELLIPSOID (KM) 
SNORM - SCHMIDT NORMALIZATION FACTORS 
Cc - GAUSS COEFFICIENTS OF MAIN GEOMAGNETIC MODEL (NT) 
CD - GAUSS COEFFICIENTS OF SECULAR GEOMAGNETIC MODEL (NT/YR) 
TC - TIME ADJUSTED GEOMAGNETIC GAUSS COEFFICIENTS (NT) 
OTIME - TIME ON PREVIOUS CALL TO GEOMAG (YRS) 
OALT - GEODETIC ALTITUDE ON PREVIOUS CALL TO GEOMAG (YRS) 
OLAT - GEODETIC LATITUDE ON PREVIOUS CALL TO GEOMAG (DEG.) 
OLON - GEODETIC LONGITUDE ON PREVIOUS CALL TO GEOMAG (DEG.) 
TIME - COMPUTATION TIME (YRS) (INPUT) 
(E.G., 1 JULY 1985 = 1985.5¢* 

ALT - GEODETIC ALTITUDE (KM (INPUT) 
GLAT - GEODETIC LATITUDE (DEG.; (INPUT) 
GLON - GEODETIC LONGITUDE (DEG.) (INPUT) 
EPOCH - BASE TIME OF GEOMAGNETIC MODEL (YRS) 

DTR - DEGREE TO RADIAN CONVERSION 
SP(M) - SINE OF (M*SPHERICAL COORD. LONGITUDE) 

CP(M) - COSINE OF (M*SPHERICAL COORD. LONGITUDE) 
ST - SINE OF (SPHERICAL COORD. LATITUDE) 
CT - COSINE OF (SPHERICAL COORD. LATITUDE) 
R - SPHERICAL COORDINATE RADIAL POSITION (KM) 
CA - COSINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE 
SA - SINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE 
BR - RADIAL COMPONENT OF GEOMAGNETIC FIELD (NT) 
BT - THETA COMPONENT OF CE *.!AGNETIC FIELD (NT) 
BP - PHI COMPONENT OF GEOMAGNETIC FIELD (NT) 
P(N,M) - ASSOCIATED LEGENDRE POL 1; NOMIALS (UNNORMALIZED) 
PP(N) - ASSOCIATED LEGENDRE POLYNOMIALS FOR M=1 (UNNORMALIZED) 
DP(N,M) - THETA DERIVATIVE OF P(N,M) (UNNORMALIZED) 
BX - NORTH GEOMAGNETIC COMPONENT (NT) 
BY - EAST GEOMAGNETIC COMPONENT (NT) 
BZ - VERTICALLY DOWN GEOMAGNETIC COMPONENT (NT) 
BH - HORIZONTAL GEOMAGNETIC COMPONENT (NT) 
DEC - GEOMAGNETIC DECLINATION (DEG.) (OUTPUT) 


EAST=POSITIVE ANGLES 
WEST=NEGATIVE ANGLES 


202 


DIP - GEOMAGNETIC INCLINATION (DEG.) (OLTPUT) 
DOWN=POSITIVE ANGLES 
UP=NEGATIVE ANGLES 
TI - GEOMAGNETIC TOTAL INTENSITY (NT) (OUTPUT) 
GV - GEOMAGNETIC GRID VARIATION (DEG.) (OUTPUT) 
REFERENCED TO GRID NORTH 


GRID NORTH REFERENCED TO 0 MERIDIAN 

OF A POLAR STEREOGRAPHIC PROJECTION 

(ARCTIC/ANTARCTIC ONLY) 
MAXDEG - MAXIMUM DEGREE OF SPHERICAL HARMONIC MODEL (INPUT) 
MOXORD ~— - MAXIMUM ORDER OF SPHERICAL HARMONIC MODEL 


ARO IO IO RR I OR IO IOI i a I OI aC OR IG IG ck kk kik ca ak ak a ak akc ak ak 


NOTE: THIS VERSION OF GEOMAG USES THE WMM-90 GEOMAGNETIC 
MODEL REFERENCED TO THE WGS-84 GRAVITY MODEL ELLIPSOID 


JCI Rie ak ak ai i ie ae keke ate eke ie tek ke ie te ie ii kk kak ki ae ki a ae ie akc ak ike ak ke ok ake ai ake ic ak ak ae ke ai ake akc akc a CC ac ake ake ake ak kc ak ake ak ake 2k 


He He Re 2c Re ae ake ate ade ake ake fe ate abe ae fe fe ake ake ae ae ate ate ake ake ale afc 2k 2h Re he ake fe ake abe abe ake fe ofc 2¥e ade ate abe ake abe abe abe ofc ofc ade oie ofc abe 24 ate abe ake ae ake ofc ake akc ake ake ake ke ai ake ake 2k 2k a ak i a 2k akc ode i 


INITIALIZATION MODULE 


ke eke ke Ne fe ke te oie ate ake ke ake ae ate fe ale ke ake ae ate ae eke ake fe ae ae fe 2Fe ote te ake oie Fe ke ake ke ie ate ake te koe ke ie aie kek ake ie ae ie ie aie ie ate ake ake ie ode te kkk ie ie ke ae a ie ee 


SUBROUTINE GEOMAG(MAXDEG) 


NA agaaaanaaNNaNnaaanananaanaaanaaananaaanaaaanaaNn 


DIMENSION C(0:12,0:12),CD(0:12,0:12),TC(0:12,0:12) 
DIMENSION P(0:12,0:12),DP(0:12,0:12),SNORM(0:12,0:12) 
DIMENSION SP(0:12),CP(0:12),FN(0:12),FM(0:12),PP(0:12) 
REAL K(0;12,0:12) 

EQUIVALENCE (SNORM,P) 


DATA EPOCH/1990.0/ 


AA AN 
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INITIALIZE CONSTANTS 


IF (MAXDEG .GT. 12) MAXDEG=12 
MAXORD=MAXDEG 
PI=3.14159265359 
DTR=PI/180.0 
SP(0)=0. 

CP(0)=1. 

P(0,0)=1. 

PP(O)=1. 

DP(0,0)=0. 
A=6378.137 
B=6356.7523142 
RE=6371.2 

A2=A**2 

B2=B**2 

C2=A2-B2 

A4=A2**2 

B4=B2**2 

C4=A4-B4 


CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO 
UNNORMALIZED 


SNORM(0,0)=1. 

DO 20 N=1,MAXORD 

ee ,0)=SNORM(N-1,0)*FLOAT(2*N-1)/FLOAT(N) 
J= 

DO 10 M=0,N 
K(N,M)=FLOAT((N-1)**2-M**2)/FLOAT((2*N-1)*(2*N-3)) 
IF (M .GT. 0) THEN 
FLNMJ=FLOAT((N-M+1)*J)/FLOAT(N+M) 
SNORMU(N,M)=SNORM(N,M-1)*SQRT(FLNMJ) 

J=1 

C(M-1,N)=SNORM(N,M)*C(M-1,N) 
CD(M-1,N)=SNORM(N,M)*CD(M-1,N) 

ENDIF 

C(N,M)=SNORM(N,M)*C(N,M) 
CD(N,M)=SNORM(N,M)*CD(N,M) 

CONTINUE 

FN(N)=FLOAT(N+1) 

FM(N)=FLOAT(N) 

CONTINUE 

K(1,1)=0. 


OTIME=- 1000. 
OALT=- 1000. 
OLAT=- 1000. 
OLON=- 1000. 


RETURN 


FOR OK ROR GOI AIO IO II GI I OK ICI I CG a IG aC ak a a ka kak ak a ok ki a a a ake ok 


PROCESSING MODULE 


SO IO IG IG IG IG II I IG GIO OR I I II i a GI a GI a ak a ak lak ica kek i akc aka ak ak ake ak a aC OK a 


ENTRY GEOMG1(ALT,GLAT,GLON, TIME,DEC,DIP,TIGV) 


AA AANAANAAAAAANANA AN 


DT=TIME-EPOCH 

IF (OTIME .LT. 0. .AND. (DT .LT. 0. .OR. DT .GT. 5.)) THEN 
PRINT *,” ’” 

PRINT *, ’ WARNING - TIME EXTENDS BEYOND MODEL S5- YEAR LIFE SPAN’ 
PRINT *, ° CONTACT DMA FOR PRODUCT UPDATES: ’ 
PRINT *,”’ 

PRINT *, ’ DEFENSE MAPPING AGENCY’ 

PRINT *, ’ ATTN: Code PRS’ 

PRINT *,’ 8613 LEE HIGHWAY’ 

PRINT *,’ FAIRFAX, VA 22031-2137’ 

PRINT *, ’ (703)285-9197, AUTOVON 356-9197’ 
PRINT *,’  ” 

PRINT *, ° EPOCH = ’,EPOCH 

PRINT *,’ TIME =’, TIME 

ENDIF 


RLON=GLON*DTR 
RLAT=GLAT*DTR 
SRLON=SIN(RLON) 
SRLAT=SIN(RLAT) 
CRLON=COS(RLON) 
CRLAT=COS(RLAT) 
SRLON2=SRLON**2 
SRLAT2=SRLAT**2 
CRLON2=CRLON**2 
CRLAT2=CRLAT**2 
SP(1)=SRLON 
CP(1)=CRLON 


CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. 


IF (ALT .NE. OALT .OR. GLAT .NE. OLAT) THEN 
Q=SQRT(A2-C2*SRLAT2) 


ore) 


aA 


ore) 


ole) 
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QI=ALT#Q 
Q2=((Q1+A2)/(Q1+B2))**2 

CT =SRLAT/SQRT(Q2*CRLAT2+SRLAT2) 
ST=SQRT(1.0-CT**2) 
R2=ALT**2+2.0*Q1+(A4-C4*SRLAT2)/Q**2 
R=SQRT(R2) 
D=SORT(A2*CRLAT2+B2*SRLAT2) 
CA=(ALT+DYR 
SA=C2*CRLAT*SRLAT/(R*D) 

ENDIF 


IF (GLON .NE. OLON) THEN 

DO 40 M=2,MAXORD 
SP(M)=SP(1)*CP(M-1)+CP(1)*SP(M- 1) 
CP(M)=CP(1)*CP(M-1)-SP(1)*SP(M-1) 
CONTINUE 

ENDIF 


AOR=RE/R 
AR=AOR**2 


BR=0. 
BT=0. 
BP=0. 
BPP=0. 


DO 70 N=1,MAXORD 
AR=AR*AOR 
DO 60 M=0,N 


COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS 
AND DERIVATIVES VIA RECURSION RELATIONS 


IF (ALT .NE. OALT .OR. GLAT .NE. OLAT) THEN 
IF (N .EQ. M) THEN 
P(N,M)=ST*P(N-1,M-1) 
DP(N,M)=ST*DP(N-1,M-1)+CT*P(N-1,M-1) 
GO TO 50 

ENDIF 

IF (N .EQ. | AND. M .EQ. 0) THEN 
P(N,M)=CT*P(N-1,M) 

DP(N,M)=CT*DP(N- 1,M)-ST*P(N-1,M) 

GO TO 50 

ENDIF 

IF (N .GT. 1 AND. N .NE. M) THEN 

IF (M .GT. N-2) P(N-2,M)=0.0 

IF (M .GT. N-2) DP(N-2,M)=0.0 
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P(N,M)=CT*P(N-1,M)-K(N,M)*P(N-2,M) 
DP(N,M)=CT*DP(N- 1,M)-ST*P(N-1,M)-K(N,M)*DP(N-2,M) 
ENDIF 
ENDIF 

50 CONTINUE 


TIME ADJUST THE GAUSS COEFFICIENTS 


IF (TIME .NE. OTIME) THEN 
TC(N,M)=C(N,M)+DT*CD(N,M) 

IF (M .NE. 0) THEN 
TC(M-1,N)=C(M-1,N)+DT*CD(M-1,N) 
ENDIF 

ENDIF 


ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS 


PAR=AR*P(N,M) 
IF (M .EQ. 0) THEN 
TEMP1=TC(N,M)*CP(M) 
TEMP2=TC(N,M)*SP(M) 

ELSE 
TEMP1=TC(N,M)*CP(M)+TC(M-1,N)*SP(M) 
TEMP2=TC(N,M)*SP(M)-TC(M-1,N)*CP(M) 
ENDIF 

BT=BT-AR*TEMP1*DP(N,M) 
BP=BP+FM(M)*TEMP2*PAR 
BR=BR+FN(N)*TEMP1*PAR 


SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES 


IF (ST .EQ. 0.0 .AND. M .EQ. 1) THEN 
IF (N .EQ. 1) THEN 

PP(N)=PP(N-1) 

ELSE 
PP(N)=CT*PP(N-1)-K(N,M)*PP(N-2) 
ENDIF 

PARP=AR*PP(N) 
BPP=BPP+FM(M)*TEMP2*PARP 
ENDIF 


AA 


60 CONTINUE 
70 CONTINUE 


ore) 


IF (ST .EQ. 0.0) THEN 
BP=BPP 
ELSE 

BP=BP/ST 


ENDIF 


ole) 
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ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO 
GEODETIC COORDINATES . 


BX=-BT*CA-BR*SA 
BY=BP 
BZ=BT*SA-BR*CA 


COMPUTE DECLINATION (DEC), INCLINATION (DIP), AND 
TOTAL INTENSITY (TI) 


BH=SQRT(BX**2+BY**2) 
TI=SQRT(BH**2+BZ**2) 
DEC=ATAN2(BY,BX)/DTR 
DIP=ATAN2(BZ,BH)/DTR 


COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT 
GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC 
(I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) 


OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0 


GV=-999.0 

IF (ABS(GLAT) .GE. 55.) THEN 
IF (GLAT .GT. 0. .AND. GLON .G 
IF (GLAT .GT. 0. .AND. GLON .L 
IF (GLAT .LT. 0. .AND. GLON .G 
IF (GLAT .LT. 0. .AND. GLON .L 
IF (GV .GT. +180.) GV=GV-360. 
IF (GV .LT. -180.) GV=GV+360. 
ENDIF 


™ 


GV=DEC-GLON 
GV=DEC+ABS(GLON) 
GV=DEC+GLON 

) GV=DEC-ABS(GLON) 


me 


0.) 
0.) 
0.) 
0. 


. 
. 


wy 


. 


OTIME=TIME 
OALT=ALT 

OLAT=GLAT 
OLON=GLON 


RETURN 


END 


CGFMFLANT (NSAP Advisor) 

CGFMFPAC (NSAP Advisor) 

SGMCDEC (NSAP Advisor) 

CINCLANTFLT (N37C, NSAP Advisor) 
CINCPACFLT (02M, 337, NSAP Advisor) 
CINCUSNAVEUR (NSAP Advisor) 
COMAREASWFORSIXTHFEFLT (CTF 66) 
COMINEWARCCM (NSAP Advisor) 
COMNAVAIRLANT (NSAP Advisor) 
COMNAVAIRPAC (NSAP Advisor) 


COMNAVSEASYSCOM (Codes 56222, Library Documentation Branch 


SEA 09B31 
COMNAVSURFLANT (NSAP Advisor) 
COMNAVSURFPAC (NSAP Advisor) 
COMOPTEVEOR 


COMPACMISTESTCEN (Codes 1018, 3250, 4024, 5021) 


COMSEABASEDASWWINGSLANT (NAVOCEANO Fit. Rep.) 
COMSECONDLET (NSAP Advisor) 

COMSEVENTHEFLT (NSAP Advisor) 

COMSIXTHELT (NSAP Advisor) 

COMSUBDEVRON TWELVE (20B) 

COMSUBLANT (NSAP Advisor) 

COMSUBPAC (NSAP Advisor) 

ALL COMSURFWARDEVGRU 

COMTHIRDFLT (NSAP Advisor) 

ALL DPT NAVSCI 


FASOTRAGRULANT (Det Brunswick, Cecil Field, Jacksonville) 
FASOTRAGRUPAC (Det Agana, Barbers Point, Cubi Point, 


Moffett Field, North Island 
FCTCLANT (Code 213) 
FLEASWTRACENLANT 
FLEASWTRACENPAC 
F LEBALMISUBTRACEN 
FLEMINEWARTRACEN 
FLENUMOCEANCEN 
NAVAIRCEVCEN (Code 8131) 
NAVAIRTESTCEN (Technical Information Dept.) 
NAVAVIONICCEN (Technical Library) 


NAVCOASTSYSCEN (Technical Information Center - Code 6120) 


NAVEAS TOCEANCEN 

NAVELEXCEN (Technical Library - Code AL) 
ALL NAVELEXDET 

NAVOCEANCOMCEN 

ALL NAVOCEANCOMDET 

NAVOCEANO (Maury Oceanographic Library) 
NAVOCEANSYSCEN (Technical Library - Code 447) 
NAVPGSCOL 

NAVPHIBASE 

NAVPOLAROCEANCEN 

NAVSHIPWPNSYSSENGSTA (Code 5125) 
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NAVSWC (Technical Library - Cede E23) 
NAVTACSUPPACT (Technical Library) 
NAVWARCOL (Technical Library) 
NAVWESTOCEANCEN 
NAVWPNSUPPCEN (Code 016) 
NISC (Technical Library - Code 63) 
NOARL. (Codes 245, 302, 352, 370, 372; 372; 550) 
NRL (Technical Library - Code 2620) 
NUSC (Technical Library - Code 02152) 
OPTEVFOR (Technical Library) 
SUBASE 
SWFPAC (Technical Library - Code SPB161) 
SWOSCOLCOM (Technical Library) 
USNA (Nimitz Library) 
WPNSTA (Technical Library) 
ARMY 
Avionics R&D Act. 
Eng Topo Lab 
ESEIA 
USAATCA-ASO 
36 Medical DET 
AIR FORCE 
AEGWC 
Ogden AFLC 
(SPECIAL) (AFLC) 
CPUSS 
HQSD 
USAF Academy (Technical Library 
9th SRW/IND 
31st Test and Evaluation Squadron 
366th Tectical Fighter Wing 
403 Rescue and Weather Reconnaissance Wing/DOX 
815 Weather Reconnaissance Wing/DOT 
4029th Strategic Reconnaissance Training Squadron 
6514th Test Squadron 
DCA (Technical Library) 
Defense Information School 
DMACSC 
DMAHTC (Tachnical Library) 
DMATAGS 
DNA (Technical Library) 
DTIC 
National Defense University 
Analysis Technology 
AVCO Systems Vision 
BGS 
Boeing Aerospace Company 
Boeing Commercial Airplane Company 
Boeing Military Airplane Company 
Canadian Pacific Airlines 
Center for Naval Analyses (Technical Library) 
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Center for Potential Field Studies, CSM 
College Observatory 

Control Systems Technology Center 
Danish Meteorological Institute 
Defense Communications Agency 
Digital Cartographic Systems 
Eastern Airlines 

EG&G/Geometrics 

ESL 

General Dynamics 

General Electric Co. 

Geophysical Services, Inc. 
Geosource Marine 

Goddard Space Flight Center 
Hughes Aircraft 

Institute Nazionale Di Geofisica 
Intergraph Ccrporation 

KLM-Royal Dutch Airlines 

LDGO 

LITEF 

Litton 

Lockheed Corporation 

Maritime Safety Agency 

McDonnell Douglas Corporation 
McDonnel Douglas Helicopter Corp. 
NOAA/NDBC 

NOAA/NGDC 

NOAA/NOS (Codes N/CG22X2, N/CG21X4) 
Northrop Corporation 

Shell Offshore, Inc. 

Singer/Link 

SIO 

Systems & Applied Sciences Corp. 
Systems Control Technology 
Technical Studies & Analytical Corp. 
TWR Corporation 

USGS Denver 

USGS Fredericksburg 

USGS Menlo Park 

WHOL 
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